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We use a one-dimensional step model to study quantitatively the growth of step bunches on Si(lll) 
surfaces induced by a direct heating current. Parameters in the model are fixed from experimental 
measurements near 900° C under the assumption that there is local mass transport through surface 
diffusion and that step motion is limited by the attachment rate of adatoms to step edges. The 
direct heating current is treated as an external driving force acting on each adatom. Numerical 
calculations show both qualitative and quantitative agreement with experiment. A force in the 
step down direction will destabilize the uniform step train towards step bunching. The average 
size of the step bunches grows with electromigration time t as t 13 , with (3 pb 0.5, in agreement 
with experiment and with an analytical treatment of the steady states. The model is extended to 
include the effect of direct hopping of adatoms between different terraces. Monte-Carlo simulations 
of a solid-on-solid model, using physically motivated assumptions about the dynamics of surface 
diffusion and attachment at step edges, are carried out to study two dimensional features that are 
left out of the present step model and to test its validity. These simulations give much better 
agreement with experiment than previous work. We find a new step bending instability when the 
driving force is along the step edge direction. This instability causes the formation of step bunches 
and antisteps that is similar to that observed in experiment. 
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I. INTRODUCTION 

In 1989 Latyshev et al^ made the startling discovery 
that a direct heating current can induce step bunching 
on vicinal Si(lll) surfaces. When the sample is resis- 
tively heated with direct current, steps can rearrange 
into closely spaced step bunches separated by wide ter- 
races. Around 900°C, the step train is unstable towards 
step bunching when the current is in the step-down direc- 
tion, but is stable when the current direction is reversed. 
Surprisingly, as the temperature is increased to 1190°, 
the stable and unstable current directions are reversed, 
i.e., the step train is unstable with step-up current and 
stable with step-down current. There is another such 
reversal as the temperature is increased further. 

Since then the phenomenon has received a great deal 
of attention. Theoretical work has mainly concen- 
trated on two goals: understanding the microscopic 
physics underlying the instability towards step bunch- 
ing and the reversal- of the unstable current direction 
with temperaturepta and determining the mesoscopic 
evolution pfi-t-he surface morphology aa-a-xesult of the 
instabilityjjliij Recently Williams et aZEjlij carried out 
a series of measurements on Si(lll) surfaces at 900°C, 
to provide a quantitative understanding of the dynam- 
ics. By controlling the experimental system and com- 
paring with theoretical models, they were able to ex- 
tract detailed information about the mechanism and 
to determine quantitative values of relevant parame- 
ters. Although the details of the microscopic mecha- 
nisms leading to the change in the destabilizing current 
direction with varying temperature are still not fully 



understood Jala we show here that there exists a reliable 
mesoscopic theory that can provide quantitative agree- 
ment with a variety of experimental results in the tem- 
perature regime (900°C) studied by Williams et al. 

In Sees. [fj and III, we briefly review some of the ex- 
perimental and theoretical work that led to our present 
model. We focus on the case where the step motion is 
limited by the attachment rate of adatoms to the step edge 
(in contrast to being limited by the diffusion rate on ter- 
races). We also assume local mass transport by surface 
diffusion. These assumptions yield a minimal mesoscopic 
model that is consistent with all previous experimental 
results. In Sec. IV we give numerical results from this 
model using realistic parameter values and interpret and 
analyze some of the results in Sec. |v|. We briefly discuss 
in Sec. VI some effects of step permeabilitylij'EJ (direct 



adatom hops from one terrace to another), which m ight 
be important in other systems, e.g., Si(001). In Sec. VII 
we present some results of Monte-Carlo simulations of 
a microscopic solid-on-solid model, using physically mo- 
tivated assumptions about the dynamics of surface dif- 
fusion and attachment at step edges. These results are 
in qualitative -agreement with experiment, in contrast to 
previous workErO using conventional Metropolis dynam- 
ics. They also help in the understanding of additional 
2D features and instabilities that cannot be described 
by the simp le ID step model. Final remarks are given 
in Sec 
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FIG. 1. Illustration of the labeling of steps and terraces 
and kinetic coefficients. 



II. ID STEP MODEL WITH EXTENDED 
VELOCITY FUNCTIONS 

Vicinal surfaces, which are created by a miscut along a 
low-index plane below the roughening temperature, are 
most naturally and usefully described by a model of in- 
teracting steps of the same sign. Because of the inherent 
anisotropy of the underlying crystal structure, these sur- 
faces often exhibit quasi one dimensional features, thus 
making a ID step model useful and accurate. Also the 
number of steps is often conserved as the surface evolves, 
permitting further simplifications in the analysis. 

The change in the morphology of vicinal surfaces can 
be described in terms of the velocity of each step as 
long as no steps are created |Or destroyed. The classic 
Burton-Cabrera-Frank (BCF)IIa treatment assumes that 
the mass transfer is governed by a set of adatom dif- 
fusion equations on each terrace, with steps acting as 
perfect sinks and sources for adatoms (or vacancies) so 
that local equilibrium is always maintained. However, 
the original BCF picture is valid only for simple materi- 
als where adatom diffusion is the rate limiting process. 
Extensions of the BCF model can be made to include 
a finite attachment/detachment term in the boundary 
conditions. This is needed for materials like the silicon, 
where the atom exchange rate between steps and terraces 
is not fast enough to permit the adatom concentration 
near the step edge to achieve local equilibrium. 

In his important work on the instability induced by a 
direct heating current on Si(lll), Stoyanovo proposed 
such a modified BCF model, including both a finite 
adatom attachment/detachment rate at step edges and 
an adatom drift velocity (or equivalently, an external 
driving force due to the electric field) . NatoriQ extended 
the work of Stoyanov to include step repulsions. The idea 
of incorporating step interactions in a generalized BCE 
model has been further developed by Sato and Uwahafl 

To describe the exchange of atoms or vacancies be- 
tween steps and their neighboring terraces (attach- 
ment/detachment), we use a linear kinetics theory and 
write the net surface flux from step n to terrace n (see 
Fig. [j] for the labeling) as 
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where fi n is the atom chemical potential at step n, de- 
fined as the increase in free energy per atom when atoms 
attach to the step. For a ID step train with elastic and 
entropic step repulsions, this can be written ascij 
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where a 2 is the area of a single atomic cell on the surface, 
h the single step height, and w n = x n+ i~x n the width of 
terrace n. Here the parameter g is just the coefficient of 
the s 3 term in the well-known Gruber-MullinsEHl form for 
the projected free energy of vicinal surfaces with slope s. 
jt/(a;„) is the adatom chemical potential on the terrace 
adjacent to step n, approaching x n from right hand (+) 
and left hand (— ) side. Ell c oq is the equilibrium adatom 
concentration on terraces. We assume there is no asym- 
metry for the kinetic coefficient for adatom attachment 
from upper and lower terraces (k+ = k_). 

To determine the adatom chemical potential on ter- 
races // (x) , we need to know the mass transport mech- 
anism on the surface. If it is much easier for adatoms to 
hop directly across a step edge from one terrace to an- 
other one in comparison to attachment at the step edge, 
then the adatom chemical potential becomes a constant 
on all terraces. We refer to this as non-local mass trans- 
port (case A) . The other limit is when there is no signifi- 
cant hopping of adatoms over the step edge, as assumed 
in the BCF model; we call this the limit of local mass 
transport (case B). Experiments on Si(lll) show that 
the relaxation rate of a step bunch with N steps scales 
with N~ a where a. ~r , 4 .3 ± 0.5. As we have discussed 
in detail elsewhere^c 2 ! this is consistent with the local 
mass transport limit (case B), and we will assume this 
limit in most of this paper. In Sec. VI we will consider 
a more general scenario. 

Assuming local mass transport, the step velocities 
v n (t) can be determined by solving the diffusion equa- 
tion for adatoms on terraces with boundary conditions 
at step edges governed by linear kinetics. The equations 
can bfi wT-ittpn generally in an extended velocity function 
/ormoa 

Vn = f+(w i) +/_(u; n _i;/x n _i,/i n ) (3) 

We will not write down the general form for the velocity 
functions f± for the electromigration problem since it is 
very complicated and not very instructive for our pur- 
pose. A simple limit that is consistent with experiment 
will be discussed below. More, general expressions have 
been given by many authors.Ej 

In studies of surface dynamics, it is often convenient 
to consider idealized models where the kinetics is limited 
by a few slow processes on the surface, and the rates 
of other faster processes are taken to infinity. For BCF 
models, neglecting evaporation and deposition, there are 
two basic rates, the attachment /detachmenL_rate k, and 
the adatom diffusion rate Bartelt et al2B estimated 
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TABLE I. Sets of parameters which give a good fit to 
the relaxation of step bunches. We use T — 2c eq a 4 K 
(ce q a 2 = 0.2ML) to compare with previous work. The ta- 
ble is taken from Ref. 15. 



Parameter Set 


d(A) 


9(e) 


r(A 3 /s) 

3 x 10 7 


D a (A 2 /s) 


A 


100000 


0.006 


5.2 xlO 11 


B 


5000 


0.006 


4 x 10 7 


3.4 xlO 10 


c 


100 


0.03 


3 x 10 s 


5.2 xlO 9 


D 


10 


0.2 


2 x 10 9 


3.5 xlO 8 



the attachment/detachment rate from the step fluctua- 
tions of Si(lll) at 900°C under the assumption that at- 
tachment/detaciupent is the rate limiting process, while 
Pimpinclli et alEl estimated from the same data the dif- 
fusion rate under the assumption that adatom diffusion 
is the rate limiting process. It is useful to define a length 
scale d = D s /k. When d is very small, the step dynamics 
is said to be diffusion limited, and when d is very large, 
the dynamics is attachment/detachment limited. 

However, direct estimation of this ratio is difficult. 
For example, Table I lists several sets of parameters 
that give good agreement with experimentsll3 on the 
electromigration-driven relaxation of step bunches on 
Si(lll), with d ranging from 10A to 10 5 A. Physically 
d has to be finite and whether a system is diffusion lim- 
ited or attachment/detachment limited depends on the 
comparison of d with other length scales, e.g., the typical 
terrace width. Here we refer to the mathematical limit 
d — ► oo as the complete attachment/detachment limited 
model, and we call the limit d — > the complete diffusion 
limited model. 

As suggested by Table I, the relaxation experiments 
can be explained using an effective charge q of the 
adatoms that reaches a finite and physically reasonable 
value as d — ► oo. In contrast, q must tend to the un- 
physical limit oo when d — ► 0. Therefore the complete 
attachment/detachment limited model is well-defined, 
and we will use this limit to illustrate the mechanism 
for electromigration. As we show below, Eq. (||) sim- 
plifies considerably in this limit. Moreover, a value of 
d > 3000A is predicted by extrapolating the diffusion 
rate from higher temperatures using a diffusion activa- 
tion energyEj of l.leV. However, numerical solutions of 
Eq. using any of the parameter sets in Table I are 
consistent with the step bunching experiments, including 
the power law for coarsening. n 

As in the Stoyanov model, □ we assume that there is 
a force F acting on each adatom because of the electric 
field. The adatom flux on terrace n under this driving 
force is 
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With complete attachment/detachment limited kinetics, 
D s tends to infinity relative to the attachment rate k, 
or j n . Therefore the adatom chemical potential ^{x) 



has a constant gradient, F, on each terrace. In gen- 
eral, fJ^(x) is affected by the motion of the neighboring 
steps, but usually the steps move very slowly so that 
they can be treated as effectively stationary as far as the 
diffusion of adatoms is concerned. Under this quasistatic 
approximation^ at any given time, the total surface flux 
into terrace n from the two neighboring steps {jn ~3n+i) 
equals to the total amount of evaporation from this ter- 
race, which is given by c cq w n /r = w n /T e , where r is the 
average lifetime of an adatom on the terrace before it 
evaporates. With these approximations, the steu-veloci- 
ties in Eq. (^) can be written in the simple formcJ'E-a 
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III. EXPERIMENTAL DETERMINATION OF 
THE PARAMETERS 

Williams et al. devised a series of experiments to mea- 
sure the various parameters and also test the assump- 
tions about the mass transport limits discussed in the 
previous section. As mentioned earlier, on Si(lll) at 
900°C, the relaxation of step bunches is consistent with 
the local mass transport limiL-The step interaction pa- 
rameter g can be measurecoEj from the distribution 
of terrace widths and step positions at equilibrium, and 
has an estimated value around 0.015eV/A 2 . Assuming 
attachment/detachment limited kinetics for mass trans- 
port, the kinetic coefficient k can be measured indepen- 
dently from the thermal fluctuations of the steps and the 
relaxation of step bunches. Rather than using k, we de- 
fine r = 2c eq a 4 n to compare with earlier work. V gives 
the step mobility for the Brownian motion of an iso- 
lated stepE3 and is measured to be around 5 x 10 7 A 3 /s. 
This value also gives a good fit to the relaxation of step 
bunches, thus providing additional evidence supporting 
the local mass transport assumption aside from the scal- 
ing behavior mentioned before. 

The force on adatoms due to the direct heating current 
can be measured from the relaxation of the step bunches 
that occuxs after reversing the current to the stable 
direction.EEl The force acting on each adatom can be con- 
veniently described in terms of an effective charge q, with 
F = qE, where the experimental value of E = 7V/cm. 
Table I lists four sets of parameters that give good fits 
to the decay of step bunches with a direct current in 
the stabilizing direction at 900°C. As mentioned before, 
as d becomes very large, the values of q and Y reach 
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FIG. 2. Evolution of average bunch size using parameter 
set A. The solid line is a fit to t 13 with f3 = 0.50. 



limiting values, which we use in the complete attach- 
ment/detachment limited model. Other relevant param- 
eters include the average terrace width wo — 1100A,|-a|nd 
the evaporation time for one monolayer r e = 1250s.E3 



IV. NUMERICAL RESULTS AND COMPARISON 
WITH EXPERIMENTS 

When the adatom drift velocity is in the step down 
direction (F < 0), one can show from a linear stability 
analysis using the parameters in Table I that a uniform 
step train is unstable towards step bunching. The evo- 
lution of the step bunches is determined by numerically 
integrating Eq. (|), starting from a step train with small 
random deviations from uniform spacing. The system 
continues to coarsen by forming larger and larger step 
bunches. Figure |^ shows the average bunch size N ave for 
a system of 2049 steps as a function of time, using Eq. (||) 
for the complete attachment/detachment limited model. 
A bunch is defined by a number of adjacent steps with 
no terraces between them larger than wq/2. The average 
bunch size is defined by J2 n n Pn/ J2 n /°«> wnere Pn is the 
number (or density) of bunches consisting of n step. It 
can be fitted by a t 13 growth rate with (3 « 0.50. Results 
using the more complicated velocity functions in Eq. (||) 
obtained from solutions of a generalized BCF equation 
and parameter set A are almost indistinguishable on this 
scale. 

This compares very well with the STM results of 
Yang, Fu, and WilliamsEl They show at both 945°C 
and 1245°C that the growth of the facet sizes Z be- 
tween two step bunches satisfies t^ where f3 » 0.5. It 
is a good approximation to relate the average number of 
steps in a bunch iV ave to Z/wq. They observed that at 
945°C, after about 120 min of annealing time, the aver- 
age terrace width between step bunches grows to about 
3500A. In the numerical simulations, the average ter- 
race width grows to about 6800A in the same time. We 
consider this quite satisfactory agreement, given the un- 
certainties in the values of the experimental parameters 



we used. Thus, not only does the step model give the 
correct power law growth rate, it also gives good quan- 
titative agreement with experimental results at 945° C. 

Simulations using other sets of parameters in Table 
I produce slightly different results, but all agree with 
the experimental data within the errors in measured pa- 
rameters. Moreover, they all have approximately a t 1 ^ 2 
coarsening rate during the time simulated and exper- 
imentally observed. Therefore we can not determine a 
unique set of microscopic parameters accurately from the 
coarsening rate alone. 

Dobbs and Krugjlil also obtained a t x l 2 coarsening 
rate from simulations of a 2D solid-on-solid model us- 
ing Metropolis dynamics. However, they obtained the 
t 1 / 2 behavior only when there is significant lateral fluc- 
tuations of step bunches as can sometimes occur in the 
later stages of coarsening, while initially the growth law 
they observed went as i 1 / 4 . Experimentally there is no 
such transition, and we obtain this kind of coarsening 
from a ID model with stra ight steps. Moreover, as dis- 
cussed below in Sec. |VH{ , there are other unphysical 
surface features that arise from the use of Metropolis 
dynamics to describe Si(lll) and related systems, and 
we suggest there an alternative dynamical scheme that 
gives good qualitative agreement with experiment. 

As another application of the step model, we also sim- 
ulate the step bunching occurring under growth pandi- 
tions. It is well known that a Schwoebel barrierE 1 ] has 
very different effects on growth and evaporation. For 
example, if there is an additional barrier for an adatom 
to attach to a step edge from the upper terrace, a ID 
step train will be stabilized under growth and destabi- 
lized under evaporation. However, even in the absence 
of a Schwoebel asymmetry (that is, even when k + = k_) 
as assumed here, simulations of the present step model 
under growth conditions show a decrease in the bunch- 
ing rate with increasing deposition rate. Figure ^ shows 
the dependence of the average bunch size as a function of 
time for different growth conditions. It is useful to define 
R = RT/c e q as the ratio between deposition and desorp- 
tion rates. As R increases, the bunching rate decreases, 
in good agreement-with the experimental results of Yang, 
Fu, and Williams x3 This decrease in the coarsening rate 
with increasing deposition was also noted by Tersoff ct 
al.E3 in their study.. of stress induced step bunching. As 
Kandel and Weeksa argued, when the step train is trav- 
eling in a certain direction (e.g., due to deposition or 
evaporation), a step at the front end of a step bunch 
can leave the bunch and join with the step bunch in 
front of it, causing an exchange of steps between step 
bunches. As the growth rate increases, the velocities of 
these crossing steps get larger and larger, so more de- 
bunching occurs, thus reducing the coarsening rate. A 
detailed study of the effect of debunching requires a 2D 
model and is beyond the scope of this paper. 



4 



100 F 



« R=0 
+ R=10/x e 
* R=56/t, 



'" 10 r 



10 



100 

time (min) 



1000 



FIG. 3. Effects of growth conditions on the step bunching. 



FIG. 4. Typical profile of step configurations when step 
bunches are induced by a force acting on adatoms in the step 
down direction. 

V. ANALYSIS: QUASI-STEADY STATE AND 
COARSENING RATE 

In this section we will try to understand analytically 
some of the numerical results. Although it is straightfor- 
ward to use the solution of the diffusion equation to de- 
termine the velocity functions in Eq. (^) when simulating 
the step bunching numerically, it is more convenient and 
instructive to consider the simple linear velocity function 
model of Eq. (||) . Initially we also neglect any deposition 
or evaporation. 

When the surface flux j is a constant everywhere, the 
surface is in a steady state. As we will show later, such 
a steady state is possible when F is in the step down 
direction. Figure [| is snapshot of the profile of the sur- 
face during one particular simulation. Assuming that 
the step bunch between A and B is in a steady state 
with a surface atom flux j* , we have from Eqs. ( |la|) and 
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where N s is the number of steps in the bunch and u>ab 
is the distance between A and B. On the flat terrace 
between B and A 1 the flux is given by 
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For a periodic array of step bunches, fiA = fJ-A' , and the 
steady state flux is then 
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Therefore the steady state adatom current is indepen- 
dent of the size of the step bunches. This result is valid 
for a periodic array of step bunches and is important 
in determining the time scaling exponent, as it will be 
shown at the end of this section. 

The steady state profile of the step bunch can be cal- 
culated numerically. Here we go to the continuum limit, 
which is a good approximation when the number of steps 
in the bunch is large. If z(x) denotes the height of 
the surface, then the slope is z x (x) = dz/dx. Using 
the contin uum version of the Gruber-Mullins free energy 
functionalcjGj 
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which is appropriate for vicinal surfaces below the ter- 
race roughening temperature, we can write the adatom 
chemical potential asE3 
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Here z xx is the second derivative of z(x) with respect to 
x. As in the step model discussed earlier, we can drop the 
linear step energy term if no new steps are created. Note 
that in the continuum description slopes of different signs 
correspond to positive and negative steps. Here we con- 
sider step profiles with positive slopes everywhere and 
thus can set \z x \ = z x . In the attachment /detachment 
limit, the surface flux is given by 
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The above equation has the physical property that the 
adatom mobility (as the response to an external field) is 
inversely proportional to the slope z x . To calculate the 
step bunch profile z{x), we can neglect the first term in 
Eq. (12) if the terrace width inside a bunch is much 
smaller than the distance between bunches. For an iso- 
lated step bunch with j(x) = j* , using Eq. (pd|), we have 
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where zq and H are integration constants. The step 
bunch profile z(x) can be easily calculated numerically 
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by integrating Eq. (|l4|), or analytically in terms of hy- 
pergeometric functions. Equation (|14| ) was first de- 
rived by NozieresES in his study of surface dynamics 
below the roughening transition. The maximum slope 
of the step bunch with j* < is given by z™ ax = 
[-j* {H/2) 2 kT I (2ngh 2 c cq )} 1/3 . For a true steady state 
with a periodic array of step bunches, j* = j* s [Eq. (Q)]. 
We expect j* to fluctuate around this value for a system 
in a quasi steady state. H is approximately the height of 

2/3 

the step bunch for large bunches. We have z™ ax ~ N s 
since j* s is independent of N s . This can be experimen- 
tally tested by measuring the average slope of the step 
bunch as a function of the average bunch size. Note that 
the continuum limit breaks down near the edges of a step 
bunch where sharp changes in the local slope occur. 

Strictly speaking, the above analysis only holds 
for steady state profiles with complete attach- 
ment/detachment limited kinetics, but we expect it to 
be a good approximation for the quasi steady state pro- 
files that arise as the step bunches slowly coarsen with 
time. Indeed simulations and experiment agree that the 
step bunches coarsen with time as t$ with (3 « 0.5. This 
t 1 / 2 power law can be justified by a scaling argument. 
We assume that as t — > oo, there is only one character- 
istic length for the system, which scales as t° . We can 
thus write all the variables on the surface in term of the 
scaled length x/t 13 at time t. We have noted that the 
steady state flux is independent of the size of the step 
bunches [Eq. ([)])]. For a system in a quasi-steady state, 
we thus assume that the flux can be written as a function 
of the scaled length only with no extra time dependence, 



i.e. 
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In contrast, the surface profile should maintain a con- 
stant average slope and thus must have the following 
scaling form: 



z{x,t)=t Z{x/t 13 ). 
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Substituting these into the equation expressing micro- 
scopic mass conservation: 



^z(x,t) ~ -J^'(M), 
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we have (3 = 1/2 by comparing the leading exponents on 
both sides of this equation. This prediction is in good 
agreement with both the experimental and the numerical 
work, as shown in Fig. 0. 



VI. EFFECTS OF STEP PERMEABILITY 

In the previous sections we focused our study of the 
step model on vicinal Si(lll) surfaces around 900° C. 
Our basic approach can be applied more generally, 



though the limits we used above are not necessarily sat- 
isfied. However as long as the mass transport is local, 
other differences from our present model, e.g., a finite dif- 
fusion length or asymmetric step edge attachment rates 
(Schwocbel barriers), can be studied using the extended 
BCF model and the velocity function approach of Eq. (||) 
in a more or less straightforward way. 

We will not detail this work here, but instead turn 
to the conceptually interesting case of step permeability, 
which can make the extended BCF picture no longer 
valid. This js motivated in rBart by recent work by 
Tanaka et alxli and Stoyanov.EJ In the classic BCF pic- 
ture with local mass transport, the only way to achieve 
adatom transport from one terrace to another one is 
through attachment to and subsequent detachment from 
the step edge separating them. However, if there is direct 
hopping of adatoms across a step edge without incorpo- 
rating into the step edge first, this causes a coupling of 
diffusion fields on adjacent terraces that must be taken 
into account. 

In a BCF-like picture, the adatom chemical poten- 
tials can have discontinuities at step edges either because 
steps are perfect sinks, or because there is a strong dif- 
fusion barrier near the step edge. In thep, analyses of 
island flattening on Si(001), Tanaka et aZ.tZI introduced 
an adatom hopping term between adjacent terraces over 
a step edge proportional to the difference between the 
chemical potentials on the two terraces. Although this 
inter-terrace hopping could be fast compared with at- 
tachment of adatoms to the step edge, it could still be 
slow compared with diffusion on flat terraces, thus al- 
lowing discontinuities in the adatom concentration field 
at step edge positions. Here we take this limit, assuming 
that the adatom diffusion rate on flat terraces is always 
much faster than both the attachment and inter-terrace 
hopping rates. 

The effect of step permeability on the relaxation of 
step bunches due to step repulsions can be studied 
straightforwardly. In the absence of an external driv- 
ing force, the adatom chemical potential on each terrace 
is a constant denoted by (see Fig. [l] for the labeling of 
terraces) . We can write the net surface flux at the right 
hand side of step n as 



On ~ « (Pn ~ Mn) + P 04-1 - 



(18a) 



and the flux at the left hand side of step n as 

On ~ « (Mn-1 ~ Mn) + P ~ Mn)- ( 18b ) 

Assuming again the quasistatic limit where j+ = 
we can solve /x^ for any given set of fj, n from a system 
of linear equations. The analytic solution is given in 
the Appendix. It is easy to see that the p — oo and 
p = limits correspond to case A and case B dynamics 
respectively. ,— ■ 

As was mentioned in Sec. H experimentsE-3 on the 
relaxation rate of step bunches on Si(lll) near 900°C 
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show a size scaling exponent a = 4.3 ± 0.5, consistent 
with the p — limit (a = 4). In comparison, if we as- 
sume p = 2k, we obtain a — 3.6 for the bunch sizes 
used in experiments and larger p will give even smaller 
a. Therefore we conclude p < 2k and can dismiss the 
importance of step permeability for Si(lll) at 900°C. 
However p can be large for other systems, or even per- 
haps for Si(lll) at different temperatures. Tanaka et 
alH estimated p = 36k on Si(001) at 950°C. Here we 
discuss some interesting effects step permeability has on 
electromigration, which could be a way to detect any 
significant permeability if it exists. 

The adatom concentration field has a constant gradi- 
ent on each terrace when there is a driving force. We 
need to be more precise in our description of the micro- 
scopic origin of step permeability to obtain a complete 
theoretical description. Here we consider the case where 
the step permeability is proportional to the difference be- 
tween the local adatom chemical potentials immediately 
to th e left and right hand sides of the step. Equations 
(18a) and (|l8b|) then become 



in ~ K [Mn - +p[[S(x n )- MnOn)L ( 19a ) 

and 

in ~ K [/AO ~ Mn] +P[^(x^) ~ Mn^n)]- (19b) 

//(x) can be determined in much the same way as be- 
fore by assuming there is a gradient, F, in /i(x) on each 
individual terrace. In a uniform step train, the exter- 
nal force creates a local chemical potential gradient, and 
introduces a "leak" of surface flux from the permeabil- 
ity term. We now show that the "leak" will create a 
long wavelength bunching instability, in contrast to the 
pairing instability familiar from the BCF picture. 

In a linear stability analysis, the step positions are 
written as 



where 



x n (t) 



u*(0) 



a in<fr+Lo(<t>)t 



u^(0) +nw , 



(20) 



(21) 



for small perturbations from uniform configurations. 
Figure ^| plots the amplification exponent w of a uniform 
step train as a function of a dimensionless wavenumber 
(<j) = 7r corresponds to the pairing mode). The solid line 
is for p/k — 100 and the dashed line is for p = 0. The 
maximum linear instability has shifted to much longer 
wavelengths. Note that very strong, repulsive interac- 
tions could also produce such a shift .0 However, for sys- 
tems with step permeability, there is a very rapid (al- 
most linear) growth in the average size of step bunches 
in the initial stage, which then crosses over to the t 1 / 2 
behavior. In contrast, for the purely repulsive system, 
the growth rate is approximately t x l 2 at all times. These 
characteristics could be used to detect step permeability 
if it is very large. 




FIG. 5. Linear instability of a uniform step train. The lo- 
cal hopping of adatoms over step edges (step permeability), 
coupled with the step repulsions, moves the maximum insta- 
bility away from the step pairing mode (q = n) to longer 
wavelengths. The solid line is for p/k = 100 and the dashed 
line is for p = 0. Other parameters are taken from set A 
in Table I, although we don't attempt to describe Si(lll) 
realistically here. 



VII. MONTE-CARLO SIMULATIONS OF THE 
2D SOLID-ON-SOLID MODEL: MODIFIED 
ARRHENIUS DYNAMICS 

For Si(lll) at 900°C the steps are mostly straight and 
a ID model is adequate for most purposes. However, 
at higher temperatures, there exists noticeable bending 
of steps. For example, at around 1100°C when evap- 
oration is significant, 2D arrays of crossing steps form 
between step bunches. Kandel and Weeksu proposed a 
(quasi) 2D step model where the velocity of each step 
depends only on the local neighboring terrace widths in 
the direction perpendicular to the average step edge di- 
rection. This model reproduced, many features of the 
crossing arrays quite accuratelyJij Further developments 
along these lines have been reported in Refs. |2^ and ^4|. 

A full 2D step model taking account of 2D adatom 
diffusion on terraces with boundary conditions on the 
moving curved steps is very difficult to study. Also it is 
necessary to go beyond the BCF framework which ex- 
cludes the creation of new steps to explain the anti-step 
bunches reported by Latyshev et alSB Here we study 
a generalized 2D solid-on-solid (SOS) model that takes 
explicit account of a step edge barrier in the kinetics of 
adatom attachment /detachment at step edges. We be- 
lieve this is probably the simplest 2D microscopic model 
that can provide a physically reasonable description of 
both adatom diffusion and step motion in Si(lll) and 
related systems. However there is almost no hope of 
simulating the long time behavior of such a microscopic 
model using realistic parameter values. Thus, in contrast 
to the ID step model we studied above, here we con- 
centrate only on qualitative properties. Specifically, we 
consider a very large external driving force along with a 
very small average terrace width. These extreme choices 
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will permit significant step motion in the computer time 
available to us. 

The SOS model is defined on a square lattice with 
total energy 

H = J2^- h j\ ( 22 ) 

(V) 

where hi is the column height and (ij) denotes nearest- 
neighbor pairs on a square lattice. Surface diffusion is 
simulated by exchange of atoms on top of a nearest- 
neighbor pair of columns (hi — > hi — 1 and hj — ► hj + 1, 
where i,j are nearest-neighbor sites). A driving force 
is simulated by asymmetric attempt frequencies in di- 
rections along and opposite the force direction. For ex- 
ample, when the force is in x direction, we assume the 
attempt frequencies in the +x and — x directions satisfy 
the following relation: 

P+x/p-x = exp(2Fa/k B T). (23) 

Our next task is to describe how the probability for 
an adatom hopping from i site to j depends on surface 
configurations. Often in statistical physics, Metropolis 
dynamics is used, where the hopping probability depends 
on the energy difference between the final and initial 
state AEij through the relation 

Ty = min[l, exp(-AK y -/fcr)]. (24) 

This dynamics often has the virtue of fast equilibration 
in the absence of a driving force since there are no barri- 
ers for movements with no change in energy, but it does 
not usually provide a physically realistic description of 
actual dynamical psqcesses. 

Krug and DobbsBiLil have studied in detail the effects 
of an external driving force combined with Metropolis 
dynamics in a SOS model. They used these simulations 
along with a continuum model to "describe the univer- 
sal features" of the electromigration problem. However, 
the resulting surface structures have several artificial fea- 
tures that do not resemble experiments on Si(lll) sur- 
faces. For example, in their simulations, a surface in- 
stability develops regardless of the current direction and 
then there are no extended flat regions of the surface 
with Vft. = 0. Experiments on Si(lll) generally reveal 
flat terraces and individual steps coalescing into bunches 
when the current is in the unstable direction, and revers- 
ing the current direction will stabilize the uniform step 
train. 

Of course it is possible that at much later times some 
limiting features of both the experiments and simula- 
tions could be insensitive to the choice of dynamics, 
and hence universal. For example, most driven sur- 
faces eventually become "rough" ; because of transverse 
step fluctuationsE3 this probably holds true in principle 
for the experiments at sufficiently large length and time 
scales even when the current is in the nominally "stable" 



direction. We show here that with a more physically mo- 
tivated choice of dynamics, the SOS model can provide a 
qualitatively accurate description of the length and time 
scales probed by present experiments, as well as of any 
longer time "universal" features, if such exist. 

Diffusion on surfaces is usually an activated process 
with an energy barrier. A different dynamical scheme, 
Arrhenius dynamics, takes this physics into account in 
an extreme way by assuming that the energy barrier is 
simply the binding energy of the atom, independent of 
the final configuration. However, Krug et alx3 found 
that there is no morphological change for this dynamics 
under an external driving force. They showed in general 
that instabilities in a continuum model are associated 
with the dependence of the adatom mobility on the lo- 
cal slope, while instabilities in a microscopic model re- 
quire a dependence of the hopping probability on the 
final configuration. Here we obtain such a configuration 
dependence by modifying the original Arrhenius dynam- 
ics, which provides a reasonable description of activated 
processes such as surface diffusion, to include an extra 
barrier that arises from the presence of steps. 

This energy barrier is motivated by the physics of re- 
bonding and surface reconstruction that can occur near 
steps. The surface atoms near steps on Si(lll) surfaces 
usually rearrange themselves and cebond in characteris- 
tic ways to lower the step energyO To incorporate an 
additional adatom into the step usually involves the col- 
lective motion of many atoms as this rebonding is mod- 
ified. This process has a higher activation energy than 
the simple pairwise additive bond picture in the usual 
SOS model would suggest. Also, in many cases the re- 
peatable step unit, the kink, has a complex structure, 
and requires the incorporation of two adatoms to bring 
about its movement. To take account of this physics in 
our simulations in a simple way, we assign an additional 
barrier for any movement that lowers the energy, since 
all attachment events are associated with a decrease in 
energy. So that detailed balance holds in equilibrium, 
the same barrier must also be added to a movement that 
increases the energy. We call this scheme modified Ar- 
rhenius dynamics and thus assume 

_ / exp(-2en l ), Ai£y = , ■. 

x « -\ bexp(-2e ni ), A£y^0, ^> 

where b < 1 and ni is the number of horizontal bonds 
the surface atom at site i has. 

This way of introducing, an attachment barrier was 
suggested by Bartelt et alxB in their study of step fluctu- 
ations. We view modified Arrhenius dynamics as a con- 
venient but not necessarily unique microscopic scheme 
that produces the "right" boundary conditions (giving 
in particular a finite value for the kinetic coefficient n) 
in the mesoscopic step models discussed in previous sec- 
tions. Thus the dynamical behavior of mesoscopic and 
macroscopic scale features in the simulations should be 
physically meaningful. 
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FIG. 6. Snapshot of a simulation using a solid-on-solid 
model with a diffusion bias perpendicular to the step edge 
direction, after about 2.1 x 10 6 Monte-Carlo steps. Param- 
eters used here are: kT = 0.8e, p x = 0.3, p- x = 0.7, 
p y = p- y = 0.5, and the size of the system is 512 x 512 
. The dark lines are normal (up) steps and white lines are 
antisteps (down steps). 

We start the simulations with a uniform step train 
with steps orientated along the vertical (y) direction. 
The height of the surface increases along the positive 
x direction. Periodic boundary conditions are used 
along the y direction. In the x direction we require 
h(x + L x ,y) = h(x,y) + No, where No is the initial 
number of steps in the system. For a system of size 
L x x Ly 1 the initial average terrace width wq = L x /No- 
With a diffusion bias in the average step down direction 
(P+x < P —xiP+y — P—y), the system is unstable towards 
step bunching. The step bunches continue to coarsen, 
consistent with the results of previous sections. Figure g 
is a snapshot of a typical configuration after some bunch- 
ing has occurred. The dark regions are step bunches, and 
single height crossing steps are visible between them. 

The qualitative features of the simulations are very 
similar to the experiments, and also to the predictions 
of the step model. Vicinal surfaces are stable during 
the time simulated when the driving force is in the step 
up direction, and unstable towards step bunching when 
the force is in the step down direction. Crossing steps 
form when there is significant evaporation. Preliminary 
results show that the coarsening rate is consistent with 
the t 1 / 2 power law, but so far the system size and sim- 
ulation time are too small to determine the exponent 
accurately. 

In the 2D step model studied in Ref. || and [24|, the 
steps are all ascending (or descending) at a given y po- 
sition. Although there is significant step bending, steps 
cannot form overhangs since step positions x n (y) are de- 
fined as single- valued functions of y. In the SOS model, 
there is no such restriction. Indeed, we can see from 
Fig. | that some crossing steps have bent so much that 




FIG. 7. The same parameters as in Fig. |6| except that 
the driving force direction is parallel to the average step edge 
direction (p y — 0.3, p- y = 0.7, p x = p- x = 0.5). The initial 
steps are along the vertical (y) direction. 

they have created anti-steps at certain y positions, i.e., 
steps of opposite sign to the initial ones at particular 
fixed y positions. In our simulations, the temperature is 
still well below the roughening temperature of the flat 
surface, but it is not energetically forbidden to create 
new steps or overhangs, in contrast to the step models 
previously studied. 

It is interesting to compar.o.these results with the ex- 
periment by Latyshev et al£B They observed anti-step 
bunch formation taking place after step bunch forma- 
tion. The first stage of the anti-step bunch formation 
occurs through the bending of the single height cross- 
ing steps between the step bunches, creating a region of 
bunched steps of the opposite sign. Indeed, we have di- 
rectly observed step bunches created from this kind of 
step bending in our model with modified Arrhenius dy- 
namics when we applied the external force in a direction 
parallel to the initial (and average) step edge direction. 
In Fig. 0, the initial (and average) step edge direction 
is in the vertical (y) direction. The bias is in the down- 
wards (—y) direction. The dark regions are step bunches 
formed by steps bending in the opposite direction to 
those individual steps (black lines) on the terraces. As 
in the previous case, following the bias (— y) direction, 
there are regions of steps going up, and regions .of anti- 
step bunches going down. We derive elsewhereEj from 
a 2D BCF-like model a new linear instability when the 
diffusion bias is parallel to the step edge direction that 
we believe underlies the patterns seen here. 

These features obtained from simulations of the new 
SOS model are quite different, from the ripple structure 
reported by Dobbs and Krugtil using Metropolis dynam- 
ics, where there are no distinct steps and facets after the 
surface develops large structures. Here the steps and 
terraces are easily discernible. Because of our more re- 
alistic treatment of the physics of surface diffusion and 
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attachment at steps and the favorable comparison with 
experiment, we believe that modified Arrhenius dynam- 
ics provides a better description for current-induced step 
bunching on Si(lll). 



VIII. CONCLUSION 

In summary, the evolution of the structure of Si(lll) 
surfaces during electromigration at 900° can be under- 
stood quantitatively using a one-dimensional step model, 
with parameter values and the mass transport mecha- 
nism determined from experiment. Specifically, the t 1 / 2 
power law growth rate for step bunch sizes is reproduced. 
We concentrate on the case where mass transport is lim- 
ited by the rate of adatom attachment to a step edge, 
but the method can be easily generalized, as illustrated 
by our discussion of direct adatom hopping between dif- 
ferent terraces. 

The ID step model has averaged over the individual 
movement of adatoms and atomic scale fluctuations of 
the steps, thus permitting simulations of the long time 
behavior using realistic parameter values. However, at 
higher temperatures, when 2D effects such as step bend- 
ing can be f^ea, even the quasi-2D step models consid- 
ered to dateucj may not be sufficient. Moreover an exact 
BCF-like treatment of full 2D diffusion problem seems 
prohibitively difficult. To examine these issues, we car- 
ried out Monte-Carlo simulations of a 2D solid-on-solid 
model, using physically motivated assumptions about 
the dynamics of surface diffusion and attachment at step 
edges. In particular we used modified Arrhenius dynam- 
ics with an extra barrier for attachment of adatoms at 
step edges and find good qualitative agreement with ex- 
periment. A new step bending instability is seen when 
there is a force acting on adatoms along the step edge 
direction that may be related to experiments by Laty- 
shev et alx3 In general we believe that this approach 
of combining information from experiment, microscopic 
simulations, and mesoscopic step models may prove use- 
ful in a number of different problems in surface science. 
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at the following system of linear equations for the terrace 
chemical potentials for a system of N steps: 
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(A2) 



and jit* = jt/( 
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). n\ can be solved for analyti- 



cally sinca-the matrix on the LHS of the equation is 
a circulantcj matrix. The result can be expressed as 
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K m describes the correlation between the adatom chemi- 
cal potential at a given terrace with the adatom chemical 
potential m steps away. K m decays exponentially as m 
increases. It is convenient to define iV c as the number of 
steps over which K m decreases by half, i.e., 



K Nc 

When p> k, we have 



K /2. 



(A6) 



(A7) 



Since N c is the range of correlation between the chemi- 
cal potential values on different terraces, the mass trans- 
port is effectively non-local over a number of steps much 
smaller than N c , and is local over a number of steps 
much larger than N c . 



APPENDIX A: 



By requiring j~ 



j„ _i_i a nd imposing periodic 



boundary conditions in Eqs. (18a) and (|18b|), we arrive 
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